The Toroidal Iron Atmosphere of a Protoneutron Star: 

Numerical Solution 

V. S. Imshennik*, K. V. Manukovskii, and M. S. Popov 

Institute for Theoretical and Experimental Physics, 
ul. Bol'shaya Cheremushkinskaya 25, Moscow, 1 17259 Russia 
Received June 20, 2003 

Abstract — A numerical method presented by Imshennik et al. (2002) is used to solve the two- 
dimensional axisymmetric hydrodynamic problem on the formation of a toroidal atmosphere during 
the collapse of an iron stellar core and outer stellar layers. An evolutionary model from Boyes et 
al. (1999) with a total mass of 25M is used as the initial data for the distribution of thermo- 
dynamic quantities in the outer shells of a high-mass star. Our computational region includes the 
outer part of the iron core (without its central part with a mass of 1M that forms the embryo of 
a protoneutron star at the preceding stage of the collapse) and the silicon and carbon— oxygen shells 
with a total mass of (1.8— 2.5)M Q . We analyze in detail the results of three calculations in which 
the difference mesh and the location of the inner boundary of the computational region are varied. 
In the initial data, we roughly specify an angular velocity distribution that is actually justified by 
the final result — the formation of a hydrostatic equilibrium toroidal atmosphere with reasonable to- 
tal mass, M tot = (0.117-0. 122)M Q , and total angular momentum, J tot = (0.445-0.472) x 10 50 erg s, 
for the two main calculations. We compare the numerical solution with our previous analytical so- 
lution in the form of toroidal atmospheres (Imshennik and Manukovskii 2000). This comparison in- 
dicates that they are identical if we take into account the more general and complex equation of 
state with a nonzero temperature and self-gravitation effects in the atmosphere. Our numerical cal- 
culations, first, prove the stability of toroidal atmospheres on characteristic hydrodynamic time scales 
and, second, show the possibility of sporadic fragmentation of these atmospheres even after a hydro- 
dynamic equilibrium is established. The calculations were carried out under the assumption of equa- 
torial symmetry of the problem and up to relatively long time scales (~10 s). © 2003 MAIK "Nau- 
ka/ Interperiodica" . 
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INTRODUCTION 

The theory of the spherically symmetric grav- 
itational collapse of iron stellar cores is known 
to have led to a clear separation of the collapse 
into two stages (see, e.g., Imshennik and Nady- 
ozhin 1982). At the first stage, the inner part of 
the iron core with a mass M iFe < 1M Q (Imshen- 
nik 1992) or, more precisely, with a mass of (0.6— 
0.8)M Q (Nadyozhin 1998) collapses homologically. 
The remaining outer part of the iron core with a 
mass M eFe = M Fe - M iFe = (0.4 - 1.4)M in the 
known range of iron core masses M Fe = (1.2 — 
2.0)M Q lags well behind its inner part in con- 
traction. In a rough approximation, it may even 
be assumed to maintain a hydrostatic equilibrium 
of the initial state. This is especially true for the 
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outer layers of a high-mass star. Imshennik et 
al. (2002) reduced the question about the hydro- 
dynamic behavior of the outer part of the iron core 
(the second stage of its collapse) and other outer 
layers of a collapsing star to the problem of their 
accretion from a hydrostatic equilibrium onto the 
embryo of a protoneutron star (PNS) with a mass 
of about 1M . The end of the first stage of iron 
core collapse was taken as the initial time of hy- 
drodynamic post-shock accretion; the corresponding 
problem was formulated by Brown et al. (1992). 
This simplification of the problem allowed us, first, 
to exclude the complex neutrino processes from 
consideration (i.e., to restrict our analysis to the 
approximation of ideal hydrodynamics), and, sec- 
ond, to include the hydrodynamic processes in the 
remaining onion stellar structure that surrounds the 
iron core of a high-mass star with a total mass 
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Mm.s > 1OM . In this paper, we continue to nu- 
merically analyze this post-shock accretion, but 
now we take into account the rotation effects by 
solving the two-dimensional axisymmetric prob- 
lem. 

The main result of the gravitational collapse of 
a rapidly rotating iron stellar core can be a binary 
system of neutron stars (Imshennik 1992) and its 
surrounding iron gaseous medium that did not en- 
ter into the composition of the rotating PNS at the 
second stage of the collapse. In this case, the ro- 
tating PNS is generally unstable against fragmen- 
tation and turns into a binary system of neutron 
stars (Aksenov et al. 1995). Based on a quasi-one- 
dimensional model, Imshennik and Nadyozhin (1992) 
showed that, with a sufficiently rapid initial rotation 
of an Mp e = 2Mq iron core and a full allowance for 
the neutrino processes, an iron atmosphere with a 
mass of (0.1— O.2)M distributed almost uniformly 
throughout the iron core with an initial outer ra- 
dius R Fe = 4.38 x 10 8 cm, i.e., with a mean density 
p eFe ^ (5.66 x 10 5 -1.13 x 10 6 ) g cm" 3 , is formed. 
The presumed long existence of this binary system 
of neutron stars, possibly for several hours (Imshen- 
nik and Popov 1994), made the problem of steady 
iron gas accretion onto this binary (actually, onto its 
more massive component located almost at the stellar 
center) of relevant interest. For the one-dimensional 
accretion of a cold (with zero temperature) iron gas 
with quasi-one-dimensional rotation, Imshennik and 
Popov (2001) managed to generalize the standard 
steady-state solution by Bondi ( 1 952) due to the very 
simple form of the specific enthalpy for the degenerate 
electron gas of an iron atmosphere with the com- 
plete ionization of iron atoms and arbitrary electron 
relativity. Application of this solution to the physical 
conditions near PNSs (with masses (1.4— 1.8)M ) 
showed that the accretion rate was very high, (20— 
3O)M s _1 at the above typical densities p eFe and a 
neutron star mass of 1.8M . On the other hand, the 

1 Imshennik et al. (2002) showed that the previously assumed 
influence of the decrease in the gravitational mass of a PNS 
due to intense neutrino energy losses during collapse was 
negligible: it does not lead to the emersion of the outer stellar 
layers in a wide range of total masses, (12— 25)M Q . Un- 
fortunately, another effect, the imitation of nucleon bubbles 
embedded in an initial spherically symmetric configuration 
under equilibrium conditions, was considered with the er- 
roneous overestimation of their total energy attributable to 
the unacceptably abrupt changes in initial specific internal 
energy at the outer boundary of the nucleon bubbles. After 
the rectification of these unfortunate inaccuracies (a several- 
fold decrease in initial specific internal energy), the emersion 
effect of the outer layers also vanished, and it was generally 
interpreted as a small correction to the adopted equation of 
state because of the energy release of free-nucleon recombi- 
nation into iron. 



steady-state accretion solution obtained does not fit 
into the iron core; besides, the time it takes for this 
solution to be established is too long. In short, it is 
not applicable to our problem on the second stage 
of the collapse of the outer iron core (Imshennik and 
Popov 2001 ). The quasi-one-dimensional rotation ef- 
fects taken into account in the above paper proved 
to be negligible because of the severe constraints 
on the specific angular momentum (at infinity!) 
specified in addition to the density p^. At the same 
time, based on ideal two-dimensional hydrodynam- 
ics, we were able to analytically construct a steady- 
state solution with the same equation of state for 
a cold iron gas in the form of toroidal atmospheres 
with a given arbitrary rotation law (Imshennik and 
Manukovskii 2000). 

A natural continuation of the above studies of 
analytical solutions seems to be the application of our 
numerical method of solution (Imshennik et al. 2002) 
to the problem of postshock accretion near neutron 
stars, using a complete equation of state with tem- 
perature effects (arbitrary electron degeneracy, the 
appearance of positrons, the presence of a nuclide gas 
and blackbody radiation), and a consistent allowance 
for the gravitational interaction, including the two- 
dimensional self-gravitation of the accreted matter. 
As was done previously (Imshennik et al. 2002), the 
hydrostatic equilibrium solutions obtained in evolu- 
tionary calculations (Boyes et al. 1999) and modified 
to incorporate the quasi-one-dimensional rotation ef- 
fects can also be used as the initial time and the 
initial conditions that describe the statement of the 
postshock accretion problem (see above). Note that 
the accuracy of the hydrostatic equilibrium of these 
outer layers, including the action of centrifugal forces 
is not of such fundamental importance as it was in 
Imshennik et al. (2002), because it essentially does 
not matter precisely what time is taken as the initial 
one in the unsteady-state hydrodynamic simulation 
of the outer layers of a high-mass star under consid- 
eration. 

As we show in this paper, the hydrodynamic 
postshock accretion actually gives rise to toroidal 
structures at the locations of the outer layers of an 
iron stellar core, in agreement with the analytical 
solutions from Imshennik and Manukovskii (2000). 
Thus, we clearly prove the hydrodynamic stability 
of such toroidal iron atmospheres against arbitrarily 
small two-dimensional perturbations. As regards the 
steady one-dimensional accretion, it does not take 
place during the postshock accretion: instead, signif- 
icant accretion onto the embryo of a PNS (with an 
initial mass of approximately 1M ) takes place in the 
near-axis region for part of the matter with the suffi- 
ciently low specific angular momenta specified in the 
initial conditions (with a total mass of about O.8M ). 
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This accretion is not only non-one-dimensional but 
also unsteady, with the ejection of whole fragments of 
matter from the toroidal atmosphere, a curious effect 
that was found in these calculations. 



FORMULATION OF THE PROBLEM 

The System of Equations, the Equation of State, 
and the Numerical Method of Solution 

In most cases, the hydrodynamic behavior of the 
envelope of a high-mass star can be described by 
the system of equations of ideal hydrodynamics. In 
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where E = e + — is the sum of the specific internal 

(e) and kinetic energies. The acceleration of gravity 
is S = Sns + gatrrp i- e -> ^ is the sum of two com- 
ponents. The first component is attributable to the 
gravitational field that is produced by the PNS located 
exactly at the coordinate origin: 

GM ns \ 

■— —,0,0). (6) 

The second component is attributable to the intrinsic 
gravitational field of the iron atmosphere. The accel- 
eration of this field is defined by the standard equation 



g atm = — V<3?, and the potential satisfies the Poisson 
equation 

A$ = 4vrG> (7) 

An efficient algorithm convenient for use in the finite- 
difference scheme of integrating the hydrodynamic 
equations on a stationary mesh in spherical coordi- 
nates (Aksenov 1999) is used to solve the Poisson 
equation (7) and determine the gravitational accel- 
eration g atm . Thus, in the problem under consid- 
eration, we abandon the Roche approximation that 
was used previously for simplicity (Imshennik and 
Manukovskii 2000; Imshennik and Popov 2001) and 
take into account the atmospheric self- gravitation 
effects. Note that, in general, they are of minor im- 
portance. 

In this paper, as in our previous paper (Imshennik 
et at. 2002), we use the equation of state for the mat- 
ter treated as a mixture of a perfect Boltzmann gas of 
nuclei with a perfect Fermi— Dirac electron— positron 
gas and blackbody radiation. In a wide temperature 
range, the matter is assumed to be a mixture of a 
baryonic component that includes free nucleons (n, 
p) and helium (|He) and iron (I^F 6 ) nuclides, and a 
Fermi— Dirac electron— positron gas with blackbody 
radiation. The equation of state obeys the nuclear 
statistical equilibrium (NSE) conditions with a con- 
stant ratio of the mass fractions of neutrons and pro- 
tons, including those bound in the helium and iron 
nuclides, 9 = 30/26 (Imshennik et at. 2002). As 
previously, we use a tabulated equation of state and 
precompute the thermodynamic functions with a high 
accuracy. This method allows complex calculations 
at each step of the solution of system (1)— (5) to be 
avoided. 

The numerical solution of our problem is based on 
the popular and universal PPM method. This method 
uses the Eulerian finite-difference scheme (Colella 
and Woodward 1984) and is a modification of Go- 
dunov's method (Godunov et at. 1976). The salient 
features of the method used were given previously 
(Imshennik et at. 2002). In particular, the real equa- 
tion of state in the numerical method is locally simu- 
lated by the so-called binomial equation of state 

P= [(j - l)e + c 2 }p - Po c 2 ; (8) 

the constants 7, Cq, and p are determined by the 
pressure and its derivative at a given point (p, e): 

't) - £ (f0 

°P J e P\ de J 

8P\ .. for 

dp 
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Using the binomial approximation for the equation 
of state considerably simplifies the solution of the 
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problem on discontinuity breakup that underlies our 
numerical method. 

Below, we also give the system of units used in our 
calculations. The scales of the physical quantities in 
this system are 

[r] = R , [V r ] = [V e ] = [V v ] = (GM /R ) 1/2 , 

(10) 

[p\ = M /(AttRI), [t\ = Rf/iGMo) 1 / 2 , 
[P] = GM 2 /(4ttR 4 ), [E] = GM /R , 
[T] = (GM 2 /(4vr^ ar)) i/4 ) 

where Rq and Mq are some length and mass scales, 
and a r is the radiation density constant. In units (10), 
the hydrodynamic equations with Newtonian gravi- 
tation contain no dimensionless parameters. It would 
be natural to use the following characteristic values 
from our problem as Rq and Mq: Rq = 10 8 cm and 
Mq = 10 32 g. The numerical values for the scales of 
the physical quantities from (10) are then 

[r] = 10 8 cm, [V r ] = [V e ] (11) 
= [V v ] = 2.583 x 10 8 cm s -1 , 
[p] = 7.958 x 10 6 gem" 3 , [t] = 3.971 x 10" 1 s, 
[P] = 5.3 10 x 10 23 erg cm" 3 , [E] = 6.674 x 10 16 erg, 
[T] = 2.894 x 10 9 K. 

Initial Conditions 

As we noted in the Introduction, the formation 
time of the PNS embryo with the mass (about 1M Q ) 
characteristic of our problem, when all of the outer 
layers may be assumed to be in hydrostatic equi- 
librium, may arbitrarily be taken as the initial time 
(t = 0). We used the distributions of thermodynamic 
quantities obtained in the studies of the evolution of 
high-mass stars (Boyes et al. 1999) as the initial 
data for constructing the initial conditions of our 
problem. From the above paper, we took the den- 
sity and temperature profiles, which, in turn, serve 
to determine the initial pressure and internal energy 
distributions by using the adopted equation of state. 
Of course, the equation of state used in our calcula- 
tions differs slightly from the equation of state used 
in obtaining these profiles. The main difference is that 
we took into account a large number of nuclides of 
chemical elements in our numerical calculations of 
stellar evolution. However, a detailed comparison of 
the equations of state (Imshennik et al. 2002) shows 
that, quantitatively, these differences are moderately 
large, and passing to a simpler equation of state in 
our calculations has no significant effect on the radial 
profiles of thermodynamic quantities. Nevertheless, 



the initial profiles were recalculated with the new 
equation of state. In these calculations, we used the 
radial density distribution from Boyes et al. (1999) 
and reconstructed the new equatorial distributions of 
the pressure (P) and the remaining thermodynamic 
quantity e by solving the system of hydrostatic equi- 
librium equations with nonzero rotation (in the form 
of a centrifugal force): 

dP (v 2 Gm\ 

^ = p {--—)> (12) 

^ = WV; (13) 

the pressure P and the radius r at the inner bound- 
ary were also taken from Boyes et al. (1999) for 
m ps IMq. Thus, strictly speaking, the initial spher- 
ically symmetric distribution of matter was one of a 
hydrostatic equilibrium only in the equatorial plane. 
The sizes of the computational region were chosen 
in such a way that the inner boundary was at the 
radius (0.876 x 10 s cm) that bounded the region of 
mass IMq in the initial profile (except the test cal- 
culation, in which the inner boundary was a factor of 
about 2 closer to the center along the radius (0.401 x 
10 8 cm)). There was arbitrariness in choosing the 
outer boundary along the radius. Since the outer, 
weakly rotating presupernova layers have no funda- 
mental effect on our result, it was not necessary to 
include them in our analysis. The outer boundary of 
the computational region in all our calculations was 
specified at a radius of 10.002 x 10 8 cm. The calcula- 
tions were performed for the evolutionary model of a 
star with a mass of 25M and solar metallicity ( Z Q = 
0.02). 

As the initial differential rotation law, we took the 
analytical formula 

u = w exp (-^> ( 14 ) 

as in our previous paper (Imshennik and Manukovskii 
2000), but as a function of the spherical radius. This 
rotation law corresponds well to the results of the 
quasi-one-dimensional hydrodynamic calculations of 
collapse by Imshennik and Nadyozhin (1992), par- 
ticularly with the given specific values of loq and r . 
We emphasize that, according to the data from Boyes 
et al. (1999), the total angular momentum Jo within 
the iron core is ~10 50 erg s; i.e., it is almost equal 
to the value from the paper mentioned above. Ex- 
pression (14) was used to construct the hydrostatic 
equilibrium distribution of thermodynamic quantities 
on the equator: in Eq. (12) was calculated by 
using the formula V v = uir. Outside the equator, the 
initial angular velocity distribution was assumed to 



THE TOROIDAL IRON ATMOSPHERE 



5 



be spherically symmetric, in accordance with (14). 
In contrast, the initial distributions of the azimuthal 
velocity and the specific angular momentum j 
were, of course, not spherically symmetric, because 
they were calculated by using the formulas 

V v = Lors'm9, j = V v rsin9. (15) 

In particular, this choice of initial data allows us to 
avoid such singularities as nonzero azimuthal velocity 
and angular momentum at 9 = 0, and softens a viola- 
tion of the hydrostatic equilibrium conditions outside 
the equator. 

Boundary Conditions 

The region of the solution of our problem, or the 
computational region, has the shape of a spherical 
envelope, r min < r < r max . The choice of r min and 
r max has already been determined by the physical con- 
siderations outlined in the Subsection "Initial Con- 
ditions." Sufficient boundary conditions should be 
specified precisely at these constant boundary values 
of the Eulerian radius. The inner boundary at r = 
r m i n is essentially a transparent wall on which the 
radial gradients of all physical quantities (V, p, and 
e) are assumed to be zero. However, it should be 
immediately noted that the fact that these derivatives 
become zero does not adequately fit the physical con- 
cept of transparency, being only its simplified, but 
formally sufficient realization. At the outer boundary, 
the boundary condition simulates a vacuum outside 
the computational region: the values of all thermody- 
namic quantities (p, e, and P) are set equal to nearly 
zero at r = r max . 

Because of equatorial symmetry in addition to ax- 
ial symmetry, we may restrict our calculations to only 
one quadrant, so the angle 6 for the computational re- 
gion ranges from to tt/2. At the 9 = ir/2 boundary, 
the velocity component V$ is assumed to be zero (and 
again the derivatives of all thermodynamic quantities 
become zero). The necessary boundary condition for 

the gravitational potential, — > as r — > oo, is 

r 

automatically satisfied in the algorithm used to solve 
the Poisson equation (Aksenov 1999). 

RESULTS OF THE NUMERICAL SOLUTION 

Below, we describe the results of our numerical 
solution of the system of differential equations (1)— 
(5) with the above initial and boundary conditions. In 
all our post-shock accretion simulations, we used the 
same model of a high-mass star with a total mass 
of 25M (Boyes et al. 1999). In the main hydrody- 
namic calculations, the outer computational bound- 
ary was exactly at a radius of 1.000168 x 10 9 cm, and 
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Fig. 1. Total mass of the matter versus time (0 < 9 < 
7r): Mi„ is the mass inside the computational region; 
Miown is the mass that passed through the inner boundary 
r — rmin! Mj P is the mass that passed through the outer 
boundary r — r max ; M gI is the total mass at r < r mm . All 
data are for calculation [2]. 

the inner computational boundary was at a radius of 
8.764082 x 10 7 cm. The cases differed in the number 
of zones into which the computational region was 
divided. In calculation [1], the total number of zones 
was 100 in the radial direction and 30 in the direction 
of change of the polar angle; in calculation [2], these 
numbers were 150 and 45, respectively. The same 
set of initial data (see the Subsection Initial Condi- 
tions) was used as the initial presupernova state. To 
solve the Poisson equation, we expanded the integral 
representation of the gravitational potential in terms 
of Legendre polynomials up to the number Z max = 20. 

The behavior of the matter was identical in all 
our calculations. The initially equilibrium (see the 
Subsection Initial Conditions) outer part of the iron 
core and the outer layers of the presupernova included 
in the computational region begin to accrete matter. 
An atmosphere in the form of a torus or, to be more 
precise, a thick disk is formed near the central PNS 
during this accretion. As a result of the accretion, 
part of the matter passes through the transparent in- 
ner computational boundary. By the time the toroidal 
atmosphere may be considered to have been formed, 
the total mass of the matter that passed inward is 
the same, with high accuracy, in all our calculations: 
~0.93M© (Fig. 1). This mass is mainly determined 
by the choice of the initial rotation law (14) for the 
matter (see uj q and r below). Because of the vacuum- 
simulating artificial boundary condition, the matter at 
the outer (along the radius) boundary also flows out 
from the computational region away from the center. 




The total mass of the matter that passes outward is 
~O.7M (Fig. 1). 

Figures 2 and 3 show lines of equal matter density 
at the final times for calculations [1] (tf = 29.034 s) 
and [2] (tf = 9.678 s), respectively. Here (and in 
Figs. 2—7, 8a, and 9—11), we use the cylindrical 
coordinates, f = rsin# and z = rcos6. Of course, 
these coordinates are equivalent to the spherical 
coordinates used in the system of equations (1)— 
(5) and the finite-difference scheme, but are more 
convenient for the physical interpretation of the 
results of our numerical solution. We see that the 
final state of the matter has an identical structure in 
both cases. The forming PNS atmosphere is almost 
twice as extended in the equatorial direction as it is 



along the rotation axis. For this reason, it would be 
more precise to call the forming atmosphere a thick 
disk rather than a torus. The density maximum in 
both calculations lies in the immediate vicinity of the 
inner computational boundary. In this situation, the 
concern about the artificial influence of the chosen 
boundary condition on the forming atmosphere is 
quite natural. This circumstance needs the following 
additional comment. 

In an equilibrium rotating atmosphere, the specific 
angular momentum of the matter is directly propor- 
tional to the fluid particle radius. Consequently, par- 
ticles with a minimum specific angular momentum 
are located near the inner boundary of the compu- 
tational region (see below for a discussion of rela- 
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r, 10 8 cm 

Fig. 3. Same as Fig. 2 for calculation [2] (t f = 9.678 s) 



tion (16)). Therefore, if we transfer the inner boundary 
to a smaller radius, then the formed region will be filled 
with matter with a lower specific angular momentum; 
thus, the density maximum will again displace to the 
inner boundary. The fact that the angular velocity in 
our initial data varies smoothly (wo = 2 s _1 , r = 5 x 
10 s cm; see the Subsection Initial Conditions), so 
that there is always an appreciable amount of matter 
with low specific angular momenta, may also con- 
tribute to this displacement. However, the rotation 
in real high-mass stars is most likely concentrated 
mainly in the region of the iron core, at the boundary 
of which the specific angular momentum decreases 
appreciably. In addition, the concentration of mat- 
ter near the inner boundary may also be partly at- 
tributable to the following factor: In our axisymmetric 



case 



d_ 

d^ 9tp 



0, the specific angular momentum 



of each fluid particle must be conserved during its mo- 
tion. However, the Eulerian scheme used here does 
not ensure the exact fulfillment of this conservation 
law. As a result, the specific angular momentum of 
fluid particles during their predominantly meridional 
flow can decrease. To check these concerns, we car- 
ried out an auxiliary calculation [3]. 

In calculation [3], which should be compared with 
calculation [ 1 ], the inner computational boundary was 
exactly at the radius of 4.010748 x 10 7 cm, and the 
distance to the radius of 8.764082 x 10 7 cm was cov- 
ered by an additional 20 zones (the computational 
region was divided into 120 cells along the radius). 
So, in the remaining part of the computational region 
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Fig. 4. Same as Fig. 2 for the test calculation [3] (tj — 6.388 s). Here, there is no other family of lines of equal angular velocity. 



(from the radius of 8.764082 x 10 7 cm to the outer 
boundary), the division into computational cells coin- 
cided exactly with the difference mesh from calcula- 
tion [1]. The minimum specific angular momentum of 
the matter from calculation [ 1 ] determined at the final 
time was taken as the minimum admissible value. For 
matter with a specific angular momentum lower than 
the value in calculation [3], this value was taken to 
be identically equal to zero. The data of this calcula- 
tion are shown in Fig. 4. The lines of equal density 
if they are carefully compared with those in Fig. 2, 
show that the coordinates of the density maximum 
remained essentially unchanged. As we see from this 
comparison, the density maximum located in Fig. 2 
near a radius r pa 0.95 is clearly established at a radius 
r pa 0.8 in Fig. 4 for the test calculation. This circum- 
stance removes the above concerns about the decisive 
influence of the boundary condition specified at the 
inner boundary on the formation of an atmosphere. 

In addition, the results obtained in calculations [1] 
and [2] should be interpreted as follows. Let j m - m be 
the minimum specific angular momentum of the at- 
mospheric matter in a steady equilibrium state. The 
simulated steady-state rotating atmosphere is then a 
result of the post-shock accretion of the outer layers 
of a high-mass star with the bulk of their rotating 
matter at the initial time having specific angular mo- 
menta higher than the above value of j min . In cal- 
culations [1] and [2], the total mass of the matter 



remaining inside the computational region by the es- 
tablishment of a steady state is ~O.13M . As we see 
from the plots in Fig. 1, a steady-state configuration 
is formed in a time of approximately 6 s; subsequently, 
the outflow of matter from the computational region 
virtually ceases, although its weak outflow results 
from the possible violation of the law of conservation 
of specific angular momentum. The rates of artificial 
accretion and ejection due to this circumstance are so 
low that they have no appreciable effect on the final 
result of our calculations. 

Figure 5 shows lines of constant angular velocity 
for the final time (tf = 9.678 s) of calculation [2]. As 
we see from the figure, the angular velocity distribu- 
tion is cylindrical in pattern; i.e., it depends only on 
the cylindrical radius f and does not depend on the 
other cylindrical coordinate z. In this case, apprecia- 
ble deviations from this rotation law appear only in the 
low-density regions outside the thick disk concen- 
trated near the axis within 6 < ir/4. The rotation law 
duj/dz = under discussion and the barotropy of the 
equation of state P = P(p) for steady-state rotating 
self-gravitating matter are known to be equivalent 
(Tassoul 1978). The lines of constant (gray) angular 
velocity were also plotted in Figs. 2 and 3. Their shape 
also confirms the aforesaid. 

More accurate data on the angular velocity distri- 
bution for calculation [2] are presented in Fig. 6. In 
this figure, the angular velocities taken on different 
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r, 10 8 cm 

Fig. 5. Lines of equal angular velocity ui at the final time for calculation [2] (tf = 9.678 s). 



rays (their angles are measured from the rotation axis) 
are plotted against the cylindrical radius f = rs'm9. 
As we see, all points satisfactorily fall on the same 
curve that is best fitted by a power law of the form uj = 
UQ(rsm6) a rather than an exponential (Gaussian) 
law of the same argument. Thus, the best fit to the 
steady-state rotation law in a thick disk is 



w = u%(r sin 6) a = 6.0562(r sinfl) 



-1.8044 -1 

^ 1 



(16) 



where the spherical radius r is given in units of the 
problem [r] = 10 s cm from (11). Note that the ex- 
ponential (Gaussian) best fit shown in Fig. 6 had 
ojq = 13.207 s" 1 and r = 1.706 (in units of [r] = 
10 s cm), but its quality was much lower than that of 
the power-law fit (16). Thus, we can easily confirm 
the above assertion that the specific angular mo- 
mentum is at a minimum near the inner boundary 
of the computational region. Indeed, j = tor 2 sin 2 6 = 
Uq(t sin#) - 1956 oc r 1956 ; i.e., the function increases 
monotonically but very slowly with radius r. 



A remark should also be made about the identi- 
fication of the presupernova envelope layers that en- 
tered into the composition of the toroidal atmosphere 
rather than collapsing onto the PNS embryo or be- 
ing ejected outward. The law of conservation of local 
specific angular momentum allows this identification 
to be made in principle, because this law is actually 
violated in the numerical solution (see above). This 
approach clearly reveals that a rather narrow layer 
of matter is transferred into the atmosphere from 
the outer part of the initial iron core and the silicon 
shell of the presupernova. The assertion that the layer 
of matter transferred into the toroidal atmosphere is 
narrow is based on the insignificant change in local 
angular momentum obtained in our calculation ( j oc 
r - 2 ). Since the change in atmospheric radius does 
not exceed a factor of 10, the angular momentum j 
changes by a factor of ~1.5 (from j min m 1.51 x 10 17 
to jmax ~ 2.28 x 10 17 cm 2 s _1 ), while the initial an- 
gular momentum changes by a factor much larger 
than ~10. Figure 7 shows the initial distribution of 




Fig. 6. Angular velocity lu versus cylindrical radius f , 
as constructed from the data on different rays (with a 
given constant angle 0) at the final time for calculation [2] 
(tf — 9.678 s): separate points on the rays adjacent to the 
equator at the location of the thick disk. The analytical 
fits are represented by the solid (power law) and dashed 
(exponential) lines. The parameters luq, lj , and ro are 
given in the text. 



Fig. 7. Lines of equal specific angular momentum j at 
the initial time for calculation [1]. The region of matter 
whose specific angular momentum lies within the range 
from j min « 1.51 x 10 17 to j max ~ 2.28 x 10 17 cm 2 s" 1 is 
colored gray; j min and j max are the limiting specific angular 
momenta in the toroidal atmosphere at the final time tf — 
29.034 s). The dotted line represents the boundary of the 
region of constant entropy. 



angular momentum j in accordance with the speci- 
fied uj distribution (14) (for the constants ujq and ro, 
see above) in spherical radius r. In this figure, the 
region of matter whose specific angular momentum j 
lies within the range from j m - m to j max (see above) and 
which, therefore, could enter the toroidal atmosphere 
is colored gray. Note that the light areas in Fig. 7 
represent the matter that either fell through the inner 
boundary to the central PNS or was ejected through 
the outer boundary into outer stellar layers. 

Let us consider the question of how strongly 
the derived nonzero temperatures of the matter in 
a toroidal atmosphere affect the numerical solution. 
The temperature reaches its maximum near the 
density maximum of a toroidal atmosphere (p max at 
''max), ~ 3 x 10 9 K, with solely iron nuclides being 
present in the nuclear composition of the atmospheric 
matter. As can be easily verified, this temperature, 
with the adopted equation of state, leads to a sig- 
nificant temperature correction to the pressure (see 
Table 1) and specific internal energy of the matter 
(see Table 2), mainly through an increase in the 
pressure and specific internal energy of the electron- 
positron matter and through the appreciable effect of 
the blackbody radiation at these temperatures. The 
estimate of the temperature effect by Imshennik and 
Manukovskii (2000) as T ~ 2 x 10 9 K also qualita- 
tively agrees with its numeric determination in our 
calculation. Let us discuss the question of whether 



these temperature corrections are consistent with the 
requirement that the equation of state for the matter 
be barotropic. Below, we discuss calculation [1]. 
Figure 8a shows the initial spherically symmetric 
distribution of specific entropy S in cylindrical ra- 
dius f on the equator. The matter that can be in 
the composition of a toroidal atmosphere (with the 
corresponding specific angular momentum) is colored 
gray. The volume in which this matter is contained 
has the shape of a hollow torus (see Fig. 7). Within 
this region, the specific entropy takes on two almost 
constant values: S\ « 2.4 x 10 8 erg g _1 Kr 1 up to 
the radius r « 3.2 x 10 8 cm (see Fig. 8a) at which 
the specific entropy abruptly changes (this radius 
is marked by the dotted line in Fig. 7) and S2 ~ 
4.4 x 10 s erg g _1 K _1 at larger radii. In Fig. 8b, the 
adiabatic curves with specific entropies Si and S2 are 
plotted in the (p, T) plane (the variables p and T are 
on a logarithmic scale). As we see from the figure, in 
the density range concerned (10 5 — 2 x 10 7 g cm -3 ), 
the simple Poisson laws with almost constant and 
similar adiabatic indices 71 « 1.346 and 72 ~ 1.307 
(see Fig. 8b) hold for S\ and S2, respectively. Thus, 
particles of matter with two different specific entropies 
(S\ and S2) could in principle enter a toroidal at- 
mosphere, which would be strange in view of the 
requirement that the equation of state for the matter of 
a toroidal atmosphere be barotropic. Interestingly, the 
total masses of these two regions in Fig. 7 are close 
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Table 1. Pressures of the matter at density p = 3.9x 10 6 g cm 3 and two different temperatures (1.0 x 10 3 and 
3.3 x 10 9 K) 



T, K 


Pad 


Pp 


Pn 


Ptot 


1.0 x 10 3 


2.500 x 10~ 3 


1.891 x 10 23 


5.791 x 10 15 


1.891 x 10 23 


3.3 x 10 9 


2.904 x 10 23 


6.028 x 10 23 


1.897 x 10 22 


9.122 x 10 23 



Note: Prad is the blackbody radiation pressure; P ep is the pressure of the Fermi— Dirac electron— positron gas; P„ is the pressure of the 
perfect Boltzmann gas of nuclides; and P tot is the total pressure of the matter (Ptot = Pad + Pep + Pn). The pressures are given in units 
of erg cm~ 3 . 



Table 2. Specific internal energies of the matter at density p — 3.9 x 10 6 gem 3 and two different temperatures (1.0 x 10 3 
and3.3xl0 9 K) 



T, K 


£rad 


£ep 


£n 


£tot 


1.0 x 10 3 


1.940 x 10~ 9 


8.540 x 10 16 


2.227 x 10 9 


8.540 x 10 16 


3.3 x 10 9 


2.382 x 10 17 


3.800 x 10 17 


7.350 x 10 15 


6.254 x 10 17 



Note: £ ra d is the specific blackbody radiation energy; £ ep is the specific energy of the Fermi— Dirac electron— positron gas; £„ is the 
specific internal energy of the perfect Boltzmann gas of nuclides; and £ tot is the total specific internal energy (£ tot = £ ra d + £e P + S n ). 
The specific internal energies are given in units of ergg -1 . 



Table 3. Physical parameters of the numerical solutions for calculations [1—3] of toroidal atmospheres 



Calculation 


[1] 


[2] 


[3] 


tf, s 


29.03449 


9.678182 


6.387759 


Pmax, 10 7 gem" 3 


0.396113 


0.374559 


1.152099 


r max , 10 8 cm 


0.955023 


0.981689 


0.766656 


M t0t , Mq 


0.117120 


0.121641 


0.211936 


J tot , 10 49 cm 2 s- 1 g 


4.453701 


4.716711 


8.152932 


E\f, 10 50 erg 


0.604490 


0.635241 


1.203042 


E^, 10 50 erg 


-1.792622 


-1.786542 


-3.893123 


10 50 erg 


0.507197 


0.528444 


1.228717 


E™1 10 50 erg 


0.501609 


0.509417 


1.169783 



Note: tf is the final time of the calculation; p max is the density maximum; r max is the cylindrical radius at which the density maximum is 
located; M tot is the total mass; J tot is the total momentum; E\° { is the total internal energy; Pg°' is the total gravitational energy; 
is the total kinetic energy; P^n is the rotational kinetic energy (oc V% /2). 



and equal to Mi = O.1136M and M 2 = O.1392M , 
respectively. Remarkably, Mi w M 2 ~ M tot (see Ta- 
ble 3) in this case. The data on the specific entropy 
distribution in the toroidal atmosphere itself, which 
are plotted in Fig. 9, introduce a larger uncertainty. 
They convincingly show that only particles with spe- 
cific entropies in the narrow range S ~ (3.3—4.5) x 
10 8 erg g _1 K _1 entered the toroidal atmosphere. 
Thus, it is clear that only the second group of particles 



with the initial specific entropy S 2 constituted the 
toroidal atmosphere obtained in our calculations. At 
the same time, the first group of particles with the 
initial specific entropy S\ was accreted onto the PNS 
embryo (or may have been ejected outward through 
the outer boundary) in the post-shock accretion 
process under study. As a result, the equation of state 
for the matter of a toroidal atmosphere proved to be 
approximately barotropic, given some errors in the 



12 



IMSHENNIK et al. 



to 

EP 

u 



x 
ii 

© 

^3 




(a) 



•(?) 



r (2) 
' max 
I 



11 



r, 10 8 



cm 



10 

10 9 t 
10 7 



6 

o 

Ml 

d.10 5 



10 3 



10 



(b) 

= Si = 2.4 x 10 8 erg g _1 K _1 (y= 1 .34568) // 

no* 

o S 2 = 4.4 x 10 s erg g"> Kr 1 (y= 1.30739) 




10 7 



10 8 



10 9 



10 



10 



T,K 



Fig. 8. (a) Specific entropy S versus cylindrical radius f 
on the equator for calculation [ 1 ] at the initial time (r^ « 



1.82 x 10 s cm, r, 



(i) 



2.35 x 10 8 cm and r 



(2) 



8.28 x 



10 s cm, rS. ~ 9.15 x 10 8 cm), (b) The adiabatic curves 
of the tabulated equation of state for two specific en- 
tropies, Si « 2.4 x 10 s erg g -1 Kr 1 and S 2 « 4.4 x 
10 8 erg g _1 K _1 . The rectilinear segments are the best 
linear fits to the adiabatic curves for the density range 
10 5 -2 x 10 7 gem" 3 . 



entropy calculation, which, according to Fig. 8b, is 
quantitatively of little importance. 

The fact that our numerical solution (see above) 
is approximately isoentropic (i.e., VS = 0) (Fig. 9) 
has a direct bearing on the dynamical stability of 
toroidal atmospheres. The Viertoft— Lebowitz stabil- 
ity criterion exists for the axisymmetric motions under 
consideration (Tassoul 1978). In our case (constant 
entropy), this criterion is 



Of 
or 



(17) 



and, as can be easily verified, it is satisfied in ac- 
cordance with (16). We immediately see from Fig. 9 
that the angular velocity distribution in regions with 
a variable specific entropy (VS / 0) deviates appre- 
ciably from the cylindrically symmetric distribution 
(dcu/dz ^ 0). This property is in close agreement 
with one of the corollaries of the more general formu- 
lation of the above criterion (Tassoul 1978). 

It is well known that an analytical solution in the 
form of rotating toroidal atmosphere can be con- 
structed in the axisymmetric case in the Roche ap- 
proximation for matter with the equation of state of 
a cold iron gas and arbitrary relativity (Imshennik 
and Manukovskii 2000). In this case, the dependence 
on the cylindrical radius alone is the only constraint 
imposed on the rotation law. If we generalize the ana- 
lytical solution from the cited paper to the power-law 
rotation ( 1 6), then for our empirically derived rotation 
law, the formula for the atmospheric density distribu- 
tion in the plane that passes through the rotation axis 
is 



p(r,9) 



£T 3 { a/i + ^Po 73 (18) 



+ 



GM 



r r ) (2a + 2)M 



x ((rsinfl)( 2a+2 ) - (r sin )( 2a + 2 )) 



3/2 



where B = 7.792 x 10~ 3 cm g" 1 / 3 and M = 2.272 x 
10 17 cm 2 s~ 2 are the constants from the equation of 
state, G is the gravitational constant, Mq is the mass 
of the central object, and (ro, Qq) are the coordinates 
of the point at which the boundary density po is spec- 
ified; the dimension relation for the coefficient loq is 
[ujo] = L^rQ J . In Fig. 10, the lines of equal density for 
this analytically specified atmosphere are compared 
with the lines of equal density for the atmosphere ob- 
tained in calculation [1]. The parameter Mo = M gr = 
1.93M (see Fig. 1) was taken from our calcula- 
tion. The constants po, ro, and 6q were chosen so 
that the total atmospheric mass from the analytical 
solution was equal to the numerically calculated at- 
mospheric mass. It should be noted that this choice 
is ambiguous. However, a common property for any 
admissible set of these parameters is not only the total 
mass of the atmosphere but also the position of its 
density maximum. The latter is determined solely by 
the properties of the analytical solution and does not 
depend on these parameters. As we see from Fig. 10, 
a more compact and denser atmosphere corresponds 
to the analytical solution. This difference could be 
explained by the use of different equations of state 
and by the fact that the meridional flow of matter, 
which is clearly present in the numerical calculation, 




Fig. 9. / — Lines of equal specific entropy at the final time for calculation [ 1 ] (tf = 29.034 s); 2 — lines of equal angular velocity w 
for the same time. 



was disregarded when deriving the analytical formula. 
In Fig. 11, atmospheric density is plotted against 
radius in the equatorial plane. This figure shows the 
density profiles for both calculations, [1, 2], taken at 
different times. As we see from the figure, the atmo- 
spheric density in calculation [2] (on a finer mesh) 
is appreciably lower than the density from calcula- 
tion [1]. The reason is as follows: in calculation [2], 
at approximately the fifth second of our calculation, 
we observed the fragmentation of the atmosphere that 
had already formed by this time; as a result, part of 
the matter was ejected outward. Concurrently with 
the fragmentation, the rotation law of the atmosphere 
changed to a form independent of the cylindrical co- 
ordinate z. In calculation [1] (on a coarse mesh), we 



observed no such event; to be more precise, weak 
perturbations were noticeable near the above time of 
the calculation (5 s), but they did not eventually lead 
to its fragmentation. We carried out an additional test 
numerical calculation [4] on an even finer mesh. The 
total number of zones in calculation [4] was 200 in 
the radial direction and 60 in the direction of change 
of the polar angle. This calculation also confirmed the 
fragmentation of the atmosphere. In both calculations 
([2, 4]), the initial perturbation that arose at the inner 
boundary of a toroidal atmosphere led to intense ejec- 
tion of a small fraction of atmospheric matter in the 
axial direction and to a division of the main volume 
of the thick disk into two parts approximately equal 
in mass. As the outer atmospheric fragment receded 
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Fig. 10. Comparison of the lines of equal atmospheric 
density at the final time for calculation [1] (/) (t/ = 
29.034 s) and the lines of equal density for the analytical 
solution (2). 
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Fig. 1 1. Atmospheric density versus cylindrical radius f 
in the equatorial plane at various times for calculations [ 1 ] 
(solid lines) and [2] (dotted lines). 

from the center, part of the matter was transferred 
back onto the remnant. As a result, the mass of 
the outward ejected matter did not exceed 30% of 
the total mass of the thick disk before the onset of 
fragmentation. The fact that fragmentation in both 
calculations (calculation [2] and test calculation [4]) 
began at the same time and developed identically 
most likely suggests that the observed fragmentation 
is not numerical in origin. 

Table 3 gives the physical model parameters ob- 
tained in calculations [1—3] of toroidal atmospheres. 
The integrated quantities in the table refer to the 
entire computational region. Without any particular 
error, they may be considered to pertain to the atmo- 
sphere, because the latter has no well-defined bound- 



ary. In addition, its surrounding matter has such a 
low density that it gives no significant contribution to 
the integrated quantities. These integrated quantities 
qualitatively agree with those obtained in the analyti- 
cal solution (cf. the data from Table 1 for M = 1.8M 
from Imshennik and Manukovskii (2000)). Of course 
in this case, there are appreciable differences between 
the differential rotation laws of the atmospheres be- 
ing compared (see above in this section). Finally, 
note that the parameters in calculation [3], in which 
the inner boundary of the computational region was 
moved to the center by twice the radius, significantly 
differ. This change can be mainly explained by the 
fact that the density maximum is closer to the center 
and corresponds to a lower specific angular momen- 
tum that a substantial mass of the matter has in the 
initial state; in particular, it can be explained by the 
shorter computational time. Unfortunately, the influ- 
ence of the boundary condition itself, which simulates 
the physical transparency condition for this Eulerian 
boundary very roughly, cannot be ruled out either. 
This question probably needs further study. 

CONCLUSIONS 

The formation of a toroidal atmosphere near a pro- 
toneutron star obtained by the hydrodynamic method 
may be considered to be the main result of our study. 
The parameters of this atmosphere with a specified 
initial rotation law in the outer layers of a high-mass 
presupernova with a structure from the evolution 
calculations by Boyes et al. (1999) were found to 
be the following: M tot = (0. 117-0. 122)M Q , J tot = 
(0.445-0.472) x 10 50 erg s, p max = (0.375-0.396) x 
10 7 g, and r max = (0.955-0.982) x 10 8 cm (cal- 
culations [1, 2]). The ranges of these quantities 
characterize the accuracy of our numerical model 
for the axisymmetric two-dimensional hydrodynamic 
problem on post-shock accretion onto the PNS 
embryo. The atmospheric parameters for the auxiliary 
calculation [3] differ markedly for the reasons given 
above (see the Section "Results of the numerical 
solution"): M tot = O.21M , J tot = 0.82 x 10 50 ergs, 
p max = 1.2 x 10 7 gem" 3 , and r max = 0.77 x 10 8 cm. 
Of course, this set of numerical results is attributable 
to the rotation law specified in the initial conditions. 
Moreover, this law provided a reasonable PNS mass 
of ~1.8M at the end of our calculation because 
of the outward ejection of a relatively large mass of 
~ O.8M and the formation of a toroidal atmosphere 
with a mass of ~O.2M (the initial mass was equal 
to the sum of ~1.OM in the PNS embryo and 
~1.8M in the outer presupernova layers up to the 
adopted radius of ~10 9 cm). This mass distribution 
between the PNS and the atmosphere was previously 
obtained in the quasi-one-dimensional calculations 
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of the collapse of an iron stellar core with a mass 
M = 2Mq by Imshennik and Nadyozhin (1992). This 
is not surprising, given the similarity between the 
rotation laws with the total angular momentum Jo ~ 
10 50 erg s specified in the initial conditions of the two 
problems (for the problem within the iron core of mass 
M Fe ~ 1.6M Q discussed here). It seems that the 
chosen rotation law has yet to be justified by gradually 
developing evolution calculations for rotating high- 
mass stars until the onset of their collapse, which is 
undoubtedly a very complex problem. The barotropic 
properties of atmospheres and the z independence of 
their angular velocity uj established in our numerical 
calculations are sufficient conditions for the validity 
of the Lichtenstein theorem (see, e.g., Tassoul 1982) 
that states the existence of an equatorial symmetry 
plane, a property that is obvious at first glance but far 
from trivial in reality. Thus, it is quite justifiable to use 
equatorial symmetry to simplify the numerical solu- 
tion (see the Subsection Boundary Conditions). 

It should also be emphasized that the numeri- 
cally calculated specific parameters of toroidal atmo- 
spheres are largely determined by the form of the 
chosen initial rotation law for the inner presuperno- 
va layers (the outer part of the iron core and the 
silicon shell). In particular, they determine the fact 
that the matter that was initially in the region of the 
silicon shell (rather than the iron core) mainly entered 
the toroidal atmosphere (see the Section "Results 
of the numerical solution"). In addition, by choos- 
ing a different initial rotation law, we could avoid 
the nonuniformity in the specific entropy distribution 
of the matter from which a toroidal atmosphere is 
formed and, as a result, avoid its fragmentation (see 
above), which may well be associated precisely with 
this distribution. Nevertheless, the very formation of 
a toroidal atmosphere during post-shock accretion is 
a universal phenomenon that depends weakly on this 
circumstance. 

The formulation of the problem on the formation of 
a toroidal atmosphere was itself largely reinforced by 
the existence of analytical solutions to the hydrostatic 
equilibrium equations for a cool iron atmosphere in 
the gravitational field of a protoneutron star (Imshen- 
nik and Manukovskii 2000). Our numerical solu- 
tion, first, shows the stability of such atmospheres 
against two-dimensional perturbations (the hydrody- 
namic sense of the relaxation method!) and, second, 
removes several restrictions in the formulation of the 
problem by Imshennik and Manukovskii (2000) due 
to the allowance for nonzero matter temperature and 



the effect of self-gravitation, which undoubtedly does 
not extend beyond minor corrections given the above 
parameters of a toroidal atmosphere. 
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